Linear Time-Invariant (LTI) Systems


By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

Table of Contents

1. Linear Time-Invariant (LTI) Systems

1.1. Systems

A discrete-time system $\mathcal{H}$ is a transformation (a rule or formula) that maps a discrete-time input signal $x$ into a discrete-time output signal $y$

1.2 Linear Systems

A system $\mathcal{H}$ is linear if it satisfies the following two properties:

1) Scaling $$\mathcal{H}\{ \alpha x \} = \alpha \mathcal{H} \quad \forall \, \alpha \in \mathbb{C}$$

2) Additivity

$$ \text{If}\,\, y_1 = \mathcal{H} \{ x_1 \} \,\, \text{and}\,\, y_2 = \mathcal{H} \{ x_2 \} \, \, \text{then } \mathcal{H} \{x_1 + x_2\} = y_1 + y_2$$

1.3 Time-Invariant Systems

A system $\mathcal{H}$ processing infinite-length signals is time-invariant (shift-invariant) if a time shift of the input signal creates a corresponding time shift in the output signal

1.4. Linear Time-Invariant (LTI) Systems

We usally consider Linear Time-Invariant (LTI) systems.

2. Convolution

Convolution is defined as the integral of the product of the two functions after one is reversed and shifted

$$y[n] = \sum \limits_{m=-\infty}^{\infty} h[n - m]\,x[m] = x[n]*h[n]$$

Output $y[n]$ came out by convolution of input $x[n]$ and system $h[n]$

  • Cross correlation and convolution
In [1]:
%%html
<center><iframe src="https://www.youtube.com/embed/Ma0YONjMZLI?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

2.1. Convolution in Matrix (= Toeplitz Matrices)

For Infinite-Length Signals

$$y[n] = x[n]*h[n] = \sum \limits_{m=-\infty}^{\infty} h[n - m]\,x[m] $$


$$ = \cdots \,h[n+2]\,x[-2] + h[n+1]\,x[-1] + h[n]\,x[0] + h[n-1]\,x[1] + h[n-2]\,x[2]\, \cdots $$


It is inner product of $h$ vectors and $x$

$$ \begin{align*} y[n] &= \begin{bmatrix} \cdots & \bracevert & \bracevert & \bracevert & \bracevert & \bracevert & \cdots\\ \cdots & h[n+2] & h[n+1] & h[n] & h[n-1] & h[n-2]\ & \cdots\\ \cdots & \bracevert & \bracevert & \bracevert & \bracevert & \bracevert & \cdots \end{bmatrix} \begin{bmatrix} \vdots\\ x[-2]\\ x[-1]\\ x[0]\\ x[1]\\ x[2]\\ \vdots \end{bmatrix}= H x \end{align*} $$

Convolution is product of matrix $H$ and $x$

$$h = h[0]$$

2.2. Impulse Response

(impulse) : $\delta[n] = \begin{cases} 1 \quad n = 0\\ 0 \quad \text{otherwise}\\ \end{cases}$

  • Output of system and delta function (impulse) is $h$. So, We call $h$ the impulse response of the system
  • From $h$, we can build matrix $H$
    • Columns/rows of $H$ are circularly shifted versions of the 0-th column/row
    • $h$ contains all the information of the LTI system
  • LTI systems are Toeplitz matrices (infinite-length signals)
    • Entries on the matrix diagonals are the same

3. Convolution Examples

Convolution = filtering in time domain

3.1. Convolution Demo

In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import IPython.display as ipd
import librosa.display

from scipy import signal

https://docs.scipy.org/doc/scipy/reference/signal.html

In [3]:
# Convolve a rectangular pulse with itself

# pulse
N = 8
n = np.arange(0, N)
x = [0,0,0,1,1,1,0,0]

# Convolve pulse with itself
y = np.convolve(x, x)

# plot
plt.figure(figsize=(10, 8))
plt.subplot(2,1,1)
plt.stem(x, markersize=4)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('pulse')

plt.subplot(2,1,2)
plt.stem(y, markersize=4)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('pulse * pulse')
plt.show()


3.2. Denoising a Piecewise Smooth Signal

  • moving average (MA) filter

  • low-pass filter in time domain

In [4]:
# piecewise smooth signal
N = 50
n = np.arange(0, N)

V = np.hstack([np.ones([1, np.int(N/2)]), -np.ones([1, np.int(N/2)])])
x = np.sin(np.pi/N*n)*V
xn = x + 0.1*np.random.randn(1, N)

# construct moving average filter impulse response of length M
M = 5
h = np.zeros([1, M])
h[0:M] = 1/M

# convolve noisy signal with impulse response
y = signal.convolve(xn, h)

# plot
plt.figure(figsize=(12,8))
plt.subplot(2,2,1)
plt.stem(x.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal')

plt.subplot(2,2,2)
plt.stem(xn.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal + noise')

plt.subplot(2,2,3)
plt.stem(h.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('impulse response')

plt.subplot(2,2,4)
plt.stem(y.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('convoluted output')
plt.show()

3.3. Edge Detection

In [5]:
# haar wavelet edge detector
M = 4
h = np.zeros([1, M])
h[0,0:2] = -1/M
h[0,2:4] = 1/M

# convolve noisy signal with impulse response
y = signal.convolve(xn, h)

# plot
plt.figure(figsize=(12,8))
plt.subplot(2,2,1)
plt.stem(x.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal')

plt.subplot(2,2,2)
plt.stem(xn.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal + noise')

plt.subplot(2,2,3)
plt.stem(h.T, 'b', 'w')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('impulse response')

plt.subplot(2,2,4)
plt.stem(y.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('convoluted output')
plt.show()
In [6]:
%%html
<center><iframe src="https://www.youtube.com/embed/I3isHB9Iy7E?rel=0"
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>
In [7]:
%%html
<center><iframe src="https://www.youtube.com/embed/ZW6pw89cuEA?rel=0" 
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>

3.4. Convolution on Audio

In [8]:
x, sr = librosa.load('./data_files/violin_origional.wav')

x = x/max(x)    # normalized

ipd.Audio('./data_files/violin_origional.wav', rate=sr) # play a wave file with sampling rate sr
Out[8]:
In [9]:
# plot wave file
plt.figure(figsize=(10, 4))
librosa.display.waveplot(x, sr=sr)
plt.xlabel('time in sec')
plt.show()
In [10]:
# impulse response in a closed room (by gunshot)

h, sr = librosa.load('./data_files/gunshot.wav')

ipd.Audio('./data_files/gunshot.wav', rate=sr) # play a wave file with sampling rate sr
Out[10]:
In [11]:
# plot wave file
plt.figure(figsize=(10, 4))
librosa.display.waveplot(h, sr=sr)
plt.xlabel('time in sec')
plt.show()
In [12]:
y = signal.convolve(x, h)
y = y/max(y)

# plot
plt.figure(figsize=(10,8))
plt.subplot(2,1,1)
librosa.display.waveplot(x, sr=sr)
plt.xlim([0, 6])
plt.title('original')

plt.subplot(2,1,2)
librosa.display.waveplot(y, sr=sr)
plt.xlim([0, 6])
plt.title('convoluted')
plt.show()
In [13]:
ipd.Audio(y, rate=sr) # play a wave file with sampling rate sr
Out[13]:

3.5. Convolution on Images

  • convolution in 2D

In [14]:
import cv2
import skimage
In [15]:
# noise image

img = cv2.imread('.\image_files\lena_sigma25.png', 0)

print(img.shape)

# plot
plt.figure(figsize = (6,6))
plt.imshow(img, cmap = 'gray')
plt.axis('off')
plt.show()

(512, 512)
In [16]:
# smoothing or noise reduction

M = np.ones([3,3])/9

img_conv = signal.convolve2d(img, M, 'valid')

# plot
plt.figure(figsize = (12,6))
plt.subplot(1,2,1)
plt.imshow(img, cmap = 'gray')
plt.title('Noisy Image')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(img_conv, cmap = 'gray')
plt.title('Smoothed Image')
plt.axis('off')
plt.show()
In [17]:
# original image

img = cv2.imread('.\image_files\lena.png', 0)

# plot
plt.figure(figsize = (6,6))
plt.imshow(img, cmap = 'gray')
plt.axis('off')
plt.show()
In [18]:
# guess what kind of image will be produced after convolution

Mv  = [[-1, 0, 1],
       [-2, 0, 2],
       [-1, 0, 1]]

Mv2 = [[1, 0, -1],
       [2, 0, -2],
       [1, 0, -1]]

img_conv = signal.convolve2d(img, Mv, 'valid')
img_conv2 = signal.convolve2d(img, Mv2, 'valid')

# plot
plt.figure(figsize = (12,6))
plt.subplot(1,2,1)
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.title('Mv')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(img_conv2, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv2))
plt.title('Mv2')
plt.axis('off')
plt.show()
In [19]:
# guess what kind of image will be produced after convolution

Mh = [[-1, -2, -1],
      [0, 0, 0],
      [1, 2, 1]]

img_conv = signal.convolve2d(img, Mh, 'valid')

# plot
plt.figure(figsize = (6,6))
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.axis('off')
plt.show()
In [20]:
M = np.array(Mv + Mh)/2

img_conv = signal.convolve2d(img, M, 'valid')

# plot
plt.figure(figsize = (6,6))
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.axis('off')
plt.show()
In [21]:
M = [[0, -1, 0],
     [-1, 4, -1],
     [0, -1, 0]]

img_conv = signal.convolve2d(img, M, 'valid')

# plot
plt.figure(figsize = (6,6))
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.axis('off')
plt.show()

4. Non-linear System: Median Filter


$$ y[n] = \text{median}\{x[n-k],\cdots,x[n+k]\}$$

  • There are nonlinear neighborhood operations that can be perfromed for the purpose of noise reduction that can do a better job of preserving edges than simple smoothing filters.

  • Median filters can do an excellent job of rejecting certain types of noise, in particular, "shot" or impulse noise (outlier in a time series) in which some individual pixels or signals have extreme values.

4.1. Removing Shot Noise in Audio

In [22]:
x, sr = librosa.load('./data_files/violin_origional.wav')

x = x/np.max(x) # normalized

ipd.Audio(x, rate=sr) # play a wave file with sampling rate sr
Out[22]:
In [23]:
# generate an audio signal with a salt and pepper noise

shot_noise = np.zeros(x.shape)
shot_noise = skimage.util.random_noise(shot_noise, mode='s&p', seed=4) 
x_noise = x + shot_noise - np.mean(shot_noise)

ipd.Audio(x_noise, rate=sr)
Out[23]:
In [24]:
# apply a linear low-pass filter

h = np.array([1/3, 1/3, 1/3])
x_avg = signal.convolve(h, x_noise)

ipd.Audio(x_avg, rate=sr)

# does not work very well
Out[24]:
In [25]:
# apply a nonlinear filter

x_median = signal.medfilt(x_noise, 7)

ipd.Audio(x_median, rate=sr)

# WOW !!!
Out[25]:
In [26]:
# plot
plt.figure(figsize=(10,16))
plt.subplot(4,1,1)
plt.plot(np.arange(0, len(x))/sr, x)
plt.ylim([-1, 1])

plt.subplot(4,1,2)
plt.plot(np.arange(0, len(x))/sr, x_noise)
plt.ylim([-1, 1])

plt.subplot(4,1,3)
plt.plot(np.arange(0, len(x)+2)/sr, x_avg)
plt.ylim([-1, 1])

plt.subplot(4,1,4)
plt.plot(np.arange(0, len(x))/sr, x_median)
plt.ylim([-1, 1])
plt.show()

4.2. Removing Shot Noise in Image

In [27]:
img = cv2.imread('.\image_files\lena.png', 0)
img_shot = skimage.util.random_noise(img, mode = 's&p')

# plot
plt.figure(figsize = (12,6))
plt.subplot(1,2,1) # original image
plt.imshow(img, cmap = 'gray', vmin = 0, vmax = np.amax(img))
plt.axis('off')

plt.subplot(1,2,2) # Degraded image with Salt & Pepper noise
plt.imshow(img_shot, cmap = 'gray', vmin = 0, vmax = np.amax(img_shot))
plt.axis('off')
plt.show()
In [28]:
# see the performance

M = np.ones([3,3])/9

img_avg = signal.convolve2d(M, img_shot, 'valid')

# plot
plt.figure(figsize = (6,6))
plt.imshow(img_avg, cmap = 'gray', vmin = 0, vmax = np.amax(img_avg))
plt.axis('off')
plt.show()
In [29]:
# Median Filter

img_median = signal.medfilt(img_shot,5)

# plot
plt.figure(figsize = (6,6))
plt.imshow(img_median, cmap = 'gray', vmin = 0, vmax = np.amax(img_median))
plt.axis('off')
plt.show()
In [30]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')